Técnicas en aprendizaje estadístico - Introduction to Statistical Learning - ISLR
Literal 4.7.1 - Ejercicio 10
a) Análisis descriptivo
El conjunto de datos Weekly contiene porcentajes de retorno semanales del índice S&P entre los años 1990 y 2010 con 1089 observaciones.
Year: es el año en el que fue obtenida la observación.Lag 1-5: Retornos porcentuales de las 5 semanas anteriores, respectivamente.Volume: Volumen de acciones intercambiadas (Número promedio de acciones diarias intercambiadas en billones).Today: Retorno porcentual de la semana.Direction: Indica si el mercado tuvo un retorno positivo o negativo en esta semana.
Encabezados y primeras 6 observaciones:
#Librerías
library(tidyverse) # Manipulación, limpieza y gráficos
library(ISLR) # Bases de datos
library(plotly) # Interactividad gráfica
library(ggthemes) # Presentación
library(GGally) # Correlación
library(yardstick) # Métricas
library(MASS) # Clasificadores
library(class) # Clasificadores
library(broom) # Resultados tidy de modelos
weekly <- as_tibble(Weekly)
head((weekly))## # A tibble: 6 x 9
## Year Lag1 Lag2 Lag3 Lag4 Lag5 Volume Today Direction
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
## 1 1990 0.816 1.57 -3.94 -0.229 -3.48 0.155 -0.27 Down
## 2 1990 -0.27 0.816 1.57 -3.94 -0.229 0.149 -2.58 Down
## 3 1990 -2.58 -0.27 0.816 1.57 -3.94 0.160 3.51 Up
## 4 1990 3.51 -2.58 -0.27 0.816 1.57 0.162 0.712 Up
## 5 1990 0.712 3.51 -2.58 -0.27 0.816 0.154 1.18 Up
## 6 1990 1.18 0.712 3.51 -2.58 -0.27 0.154 -1.37 Down
Volumen promedio de acciones diarias intercambiadas en el tiempo
Se evidencia un crecimiento constante con un alza grande en el año 2005 y una estabilización de la trayectoria alrededor del año 2007. Se añadió un ruido a la posición de cada punto para poder evidenciar la densidad, además de transparencia. Se utilizó el método mgcv::gam para ajustar la tendencia.
p <- weekly %>%
ggplot(aes(x = Year, y = Volume)) +
geom_jitter(alpha = 0.6) +
geom_smooth() +
labs(x = NULL,
y = "Volumen (Billones)") +
theme_economist()
p_box <- weekly %>%
ggplot(aes(x = as_factor(Year), y = Volume, color = Today)) +
geom_boxplot() +
theme_economist() +
theme(axis.text.x = element_text(angle = 45)) +
labs(x = NULL,
y = "Volumen (Billones)")
ggplotly(p, width = 800, height = 500)El incremento en la volatilidad del índice acompañó su crecimiento a partir del 2007.
Correlación entre las variables
No se encontraron relaciones lineales claras entre las variables. Se presentan comportamientos esperados de la variable Today con la variable categórica Direction representada en el color, donde los valores negativos están asociados a un decrecimiento.
b) ¿Podemos predecir si el índice subirá o no mediante regresión logística?
Utilicemos regresión logística para predecir si habrá un retorno positivo del mercado en esta semana. Para ello, se utilizará el retorno porcentual de las 5 últimas semanas y el volumen sin separar un conjunto de entrenamiento y de prueba (Modelo sobreajustado).
market_logistic <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume, data = weekly, family = binomial)
summary(market_logistic)##
## Call:
## glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 +
## Volume, family = binomial, data = weekly)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6949 -1.2565 0.9913 1.0849 1.4579
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.26686 0.08593 3.106 0.0019 **
## Lag1 -0.04127 0.02641 -1.563 0.1181
## Lag2 0.05844 0.02686 2.175 0.0296 *
## Lag3 -0.01606 0.02666 -0.602 0.5469
## Lag4 -0.02779 0.02646 -1.050 0.2937
## Lag5 -0.01447 0.02638 -0.549 0.5833
## Volume -0.02274 0.03690 -0.616 0.5377
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1496.2 on 1088 degrees of freedom
## Residual deviance: 1486.4 on 1082 degrees of freedom
## AIC: 1500.4
##
## Number of Fisher Scoring iterations: 4
Si realizamos una prueba de hipótesis formal para la significancia de los parámetros, encontraremos que la probabilidad de rechazar la hipótesis nula es muy grande para todas las variables predictoras excepto para el Lag2 donde al parecer hay significancia. El valor en el tiempo \(t-2\) parece tener una relación con el tiempo \(t\).
c) Matriz de confusión de la regresión logística y calidad de la predicción del modelo sobreajustado
## Up
## Down 0
## Up 1
market_logistic_class <- tibble(value = predict(market_logistic, type = "response"),
pred = fct_relevel(as_factor(case_when(value > 0.5 ~ "Up",
TRUE ~ "Down")), "Down", "Up"),
real = weekly$Direction)
market_logistic_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo de regresión logística")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.561
## 2 kap binary 0.0350
Se obtiene una tasa de clasificación correcta accuracy del 56.1% utilizando regresión logística. Es importante interpretar este resultado correctamente, ya que la predicción obtenida indica que el mercado subió 557 días y bajó 54 días, y que no se está prediciendo el valor en el tiempo \(t+1\).
Cada predicción individual se debe interpretar de la siguiente manera: “Dadas las 5 semanas anteriores y el volumen de ésta, el mercado tendrá un retorno positivo”. La mayor cantidad de errores se cometió prediciendo que el mercado iba a subir y en realidad bajó, esta situación ocurrió 430 veces.
d) Modelo de regresión logística sólo con Lag2 y conjunto de prueba
weekly_train <- weekly %>%
filter(Year <= 2008)
weekly_test <- weekly %>%
filter(Year > 2008)
lag2_logistic <- glm(Direction ~ Lag2, data = weekly_train, family = "binomial")
summary(lag2_logistic)##
## Call:
## glm(formula = Direction ~ Lag2, family = "binomial", data = weekly_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.536 -1.264 1.021 1.091 1.368
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.20326 0.06428 3.162 0.00157 **
## Lag2 0.05810 0.02870 2.024 0.04298 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1354.7 on 984 degrees of freedom
## Residual deviance: 1350.5 on 983 degrees of freedom
## AIC: 1354.5
##
## Number of Fisher Scoring iterations: 4
lag2_class <- tibble(value = predict(lag2_logistic, type = "response", newdata = weekly_test),
pred = fct_relevel(as_factor(case_when(value > 0.5 ~ "Up",
TRUE ~ "Down")), "Down", "Up"),
real = weekly_test$Direction)
lag2_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para la regresión logística con Lag2 como único predictor")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.625
## 2 kap binary 0.141
Al igual que el modelo anterior, se tienen fallas al intentar predecir las bajas del mercado, sin embargo, este se validó en un conjunto de prueba de 104 observaciones desconocidas para el modelo. El accuracy en este caso es del 62.5%. Los resultados son sorpresivamente mejores que en el modelo sobreajustado del literal b).
d) Modelo LDA (Linear Discriminant Analysis) para predecir la tendencia del mercado
## Call:
## lda(Direction ~ Lag2, data = weekly_train)
##
## Prior probabilities of groups:
## Down Up
## 0.4477157 0.5522843
##
## Group means:
## Lag2
## Down -0.03568254
## Up 0.26036581
##
## Coefficients of linear discriminants:
## LD1
## Lag2 0.4414162
lag2_lda_class <- tibble(pred = predict(lag2_lda, newdata = weekly_test)$class,
real = weekly_test$Direction)
lag2_lda_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo LDA con Lag2 como único predictor")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.625
## 2 kap binary 0.141
Se obtienen las mismas clasificaciones para el modelo LDA y el modelo de regresión logística. Como el coeficiente discriminante lineal es Lag2 = 0.4414, es de esperar que en ambas categorías la gráfica de éstos sea similar. La matriz de confusión es igual a la del modelo logístico.
f) Modelo QDA (Quadratic Discriminant Analysis) para predecir la tendencia del mercado
## Call:
## qda(Direction ~ Lag2, data = weekly_train)
##
## Prior probabilities of groups:
## Down Up
## 0.4477157 0.5522843
##
## Group means:
## Lag2
## Down -0.03568254
## Up 0.26036581
lag2_qda_class <- tibble(pred = predict(lag2_qda, newdata = weekly_test)$class,
real = weekly_test$Direction)
lag2_qda_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo QDA con Lag2 como único predictor")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.587
## 2 kap binary 0
El clasificador QDA clasificó todos los valores como “Up”, si bien es una estrategia que otorga resultados buenos en términos porcentuales con respecto al accuracy, se espera que el modelo proponga en algunos casos predicciones de baja del mercado.
g) Modelo KNN (K-Nearest Neighbors) para predecir la tendencia del mercado
set.seed(1)
lag2_knn <- knn(train = as.matrix(weekly_train$Lag2), test = as.matrix(weekly_test$Lag2), cl = weekly_train$Direction, k = 3)
lag2_knn_class <- tibble(pred = lag2_knn,
real = weekly_test$Direction)
lag2_knn_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo KNN con Lag2 como único predictor")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.548
## 2 kap binary 0.0453
El clasificador \(KNN\) tiene una tasa de clasificación correcta del 50% utilizando \(k = 1\), sin embargo, al utilizar \(k = 3\) se obtiene un 54.8%. Este clasificador no tuvo inconvenientes en proponer predicciones variadas según el valor anterior del mercado.
h) ¿Qué métodos producen los mejores resultados para predecir la tendencia del mercado utilizando Lag2?
Los métodos de regresión logística y análisis de discriminante lineal obtuvieron los mejores resultados, con un accuracy del 62.5%. Se decide favorecer al modelo de regresión logística dada la facilidad de interpretación de los coeficientes y las facilidades gráficas de comunicación que permite la curva logística.
g <- augment_columns(lag2_logistic, type.predict = "response", newdata = weekly_test) %>%
arrange(.fitted) %>%
rownames_to_column() %>%
mutate(rowname = as.integer(rowname)) %>%
inner_join(lag2_class, by = c(.fitted = "value")) %>%
ggplot(aes(x = rowname, y= .fitted, color = Direction, label = Year, label1 = pred)) +
geom_point(alpha = 0.5, shape = 1, stroke = 2) +
geom_hline(yintercept = 0.5) +
theme_light() +
xlab("Índice") +
ylab("Probabilidad de que el mercado suba")
ggplotly(g, width = 800, height = 500)i) Mejora del mejor modelo y variables adicionales de interés
Selección Step-wise con la función stepAIC sin interacciones
full_logistic <- glm(Direction ~ .-Year-Today, data = weekly_train, family = "binomial")
stepAIC(full_logistic, direction = "both", trace = FALSE)##
## Call: glm(formula = Direction ~ Lag1 + Lag2, family = "binomial", data = weekly_train)
##
## Coefficients:
## (Intercept) Lag1 Lag2
## 0.21109 -0.05421 0.05384
##
## Degrees of Freedom: 984 Total (i.e. Null); 982 Residual
## Null Deviance: 1355
## Residual Deviance: 1347 AIC: 1353
step_logistic <- glm(Direction ~ Lag1 + Lag2, data = weekly_train, family = "binomial")
step_class <- tibble(value = predict(step_logistic, type = "response", newdata = weekly_test),
pred = fct_relevel(as_factor(case_when(value > 0.5 ~ "Up",
TRUE ~ "Down")), "Down", "Up"),
real = weekly_test$Direction)
step_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para la regresión logística")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.577
## 2 kap binary 0.0350
El modelo seleccionado por Stepwise selection tiene un accuracy menor que el de un solo predictor, sin embargo, su AIC es el mismo. El AIC evalúa únicamente el ajuste, por lo que se recomendaría el modelo de una sola variable predictora. Es importante resaltar que el modelo de menor AIC podría ser mejor en un conjunto de prueba diferente.
Mejor modelo por tanteo de interacciones
Se proponen variables con interacciones y se selecciona el modelo con mayor accuracy en el conjunto de prueba, en este caso utilizando la variable Lag2 y la interacción entre Lag1, Lag 4 y Volume
final_logistic <- glm(Direction ~ Lag2 + Lag1:Lag4:Volume, data = weekly_train, family = "binomial")
summary(final_logistic)##
## Call:
## glm(formula = Direction ~ Lag2 + Lag1:Lag4:Volume, family = "binomial",
## data = weekly_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.493 -1.265 1.023 1.089 1.446
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2047966 0.0643540 3.182 0.00146 **
## Lag2 0.0553928 0.0291831 1.898 0.05768 .
## Lag1:Lag4:Volume 0.0007316 0.0014747 0.496 0.61982
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1354.7 on 984 degrees of freedom
## Residual deviance: 1350.3 on 982 degrees of freedom
## AIC: 1356.3
##
## Number of Fisher Scoring iterations: 4
final_class <- tibble(value = predict(final_logistic, type = "response", newdata = weekly_test),
pred = fct_relevel(as_factor(case_when(value > 0.5 ~ "Up",
TRUE ~ "Down")), "Down", "Up"),
real = weekly_test$Direction)
final_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para la regresión logística")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.644
## 2 kap binary 0.179
Con base en los resultados anteriores se selecciona el modelo por tanteo por su accuracy superior, sin embargo, esto puede ser un sobreajuste del conjunto de prueba; por ello, si este modelo fuera a producción se recomienda realizar validación cruzada para identificar si la interacción de variables añadida es significativa. Si no es lo es, el modelo seleccionado por Stepwise selection se considera apropiado.